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Abstract 

We analytically study the input-output properties of a neuron whose active dendritic tree, mod- 
eled as a Cayley tree of excitable elements, is subjected to Poisson stimulus. Both single-site and 
two-site mean-field approximations incorrectly predict a non-equilibrium phase transition which is 
not allowed in the model. We propose an excitable- wave mean-field approximation which shows 
good agreement with previously published simulation results [Gollo et al, PLoS Comput. Biol. 
5(6) el000402 (2009)] and accounts for finite-size effects. We also discuss the relevance of our re- 
sults to experiments in neuroscience, emphasizing the role of active dendrites in the enhancement 
of dynamic range and in gain control modulation. 
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I. INTRODUCTION 

Computational neuroscience is a growing field of research which attempts to incorporate 
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Since the 



increasingly detailed aspects of neuronal dynamics in computational models 
pioneering work of Hodgkin and Huxley (HH) jsl , which unveiled how the action potential 
in the giant squid axon could be described by ordinary differential equations governing 
the gating of ionic conductances across a membrane patch, the computational modeling of 
neuronal biophysical processes has been done at several levels, from whole neural networks 
to dendritic spines and even single ionic channel dynamics jj]. 

Rail was probably the first to extend conductance-based modeling to dendrites , start- 
ing what is nowadays a field of its own: the investigation of so-called dendritic computa- 



tion [6|. The main theoretical tool in this enterprise has been cable theory, the extension 
[via partial differential equations (PDEs)] of the HH formalism to extended systems, which 
allows one to include spatial information about dendrites such as the variation of channel 
densities along the trees, different branching patterns etc. Q]- The assumption that den- 
drites are passive elements renders cable theory linear, allowing the application of standard 
techniques from linear PDEs and yielding insightful analytical results [7|. This assumption, 
however, has been gradually revised since the first experimental evidences that dendrites 
have nonlinear properties 8|. A variety of channels with regenerative properties are now 
identified which can sustain the propagation of nonlinear pulses along the trees (called den- 
dritic spikes), whose functional role has nonetheless remained elusive 

The conditions for the generation and propagation of dendritic nonlinear excitations 
have been investigated via cable theory at the level of a dendritic branchlet. This 

has proven useful for understanding the specific role of each ionic channel in the dynamical 
properties of the nonlinear propagation, specially in comparison with experiments, which 
have mostly been restricted to the injection of current at some point in the neuron (say, 
a distal dendrite) and the measurement of the membrane potential at another point (say. 



the soma) 10|, While this limitation is justified by the difficulties of injecting currents 
and measuring membrane potentials in more than a couple of points in the same neuron, we 
must remember that neurons in vivo are subjected to a different stimulus regime, with many 
synaptic inputs arriving with a high degree of stochasticity and generating several dendritic 
spikes which may propagate and interact. 

In this more realistic and highly nonlinear scenario, cable theory, though still having the 
merit of being able to incorporate as many ionic channels as experiments reveal, becomes 
analytically untreatable. Being able to reproduce the fine-grained experimental results of a 
complex system such as a neuron does not imply that the essential aspects of its dynamics 
will be identified. Or, to put it in a renormalization group parlance, "realistic biophysical 
modeling" does not allow us to separate the relevant observables from the irrelevant ones 
that can be eliminated without significantly changing some robust property of the system. 
In fact, this has been recognized in the neuroscience literature, which has emphasized the 



need for theoretical support 12 
field of dendritic computation 



ljl9|. 



and witnessed the increase of theoretical papers in the 



In this context, we have recently attempted to understand the behavior of an active 



dendritic tree by modeling it as a large network of interacting nonlinear branchlets un- 
der spatio-temporal stochastic synaptic input and allowing for the interaction of dendritic 
spikes 20j. With a statistical physics perspective in mind, we have tried to incorporate in 
the model of each branchlet only those features that seemed most relevant, and have inves- 
tigated the resulting collective behavior. Thus each excitable branchlet was modeled as a 
simple 3-state cellular automaton, with the propagation of dendritic spikes occurring with 
probabilities which depend on direction (to account for the differences between forward- and 
backward-propagating spikes). 

This model has revealed that such a tree performs a highly nonlinear "computation", 
being able to compress several decades of input rate intensity into a single decade of out- 
put rate intensity. This signal compression property, or enhancement of dynamic range, is 
a general property of excitable media and has proven very robust against variations in the 
topology of the medium and the level of modeling, from cellular automata to compartmental 



conductance-based models 



2lN33|. Furthermore, the idea that dynamic range can be en- 



hanced in neuronal excitable media has received support from experiments in very different 
setups [34, l35|, which again suggests that the phenomenon is robust. 



Our aim here is to analytically explore the model introduced in Ref. (20| and described 
in section |Tll In section IIIII we show that the traditional cluster approximations applied to 
the system master equations fail to qualitatively reproduce the essential features observed 
in the simulations and experimental data. We propose a mean-field approximation which 
circumvents the problems faced by the traditional approach, yielding good agreement with 
simulations. We conclude in section HV] with a discussion of the consequences of our results 
for neuroscience and the perspectives for future work. 



II. MODELING AN ACTIVE DENDRITIC TREE 

The dendritic tree of an isolated neuron contains no loops and divides in two daughter 
branches at branching points. For instance. Fig. [T] (a) depicts one of Ramon y Cajal's 
drawings of a human Purkinje cell, which shows a huge ramification. Measured by the 
average number G of generations (i.e., the number of branch-doubling iterations the primary 
dendrite undergoes), the size of the dendritic trees can vary widely. One can think of an 
active dendritic tree as an excitable medium jsol , in which each site represents, for instance, 
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a branching point or a dendritic branchlet connected with two similar sites from a higher 
generation and one site from a lower generation. Correspondingly, the standard model in this 
paper is a Cayley tree with coordination number z = 3 



20|. Each site at generation g has a 
mother branch from generation g — 1 and generates two daughter branches [k = z — 1 = 2) 
at generation g + 1. The single site at g = would correspond to the primary (apical) 
dendrite which connects with the neuron soma [see Fig. [11(b)]. Naturally, the Cayley tree 
topology of our model is a crude simplification of a real tree, as attested by the differences 
between Figs. [1] (a) and[T]^b). 





FIG. 1. (Color online) Model of an active dendritic tree, (a) A famous drawing by Ramon y 
Cajal of a human Purkinje cell, (b) Excitable elements (circles) connected (bars) in a Cayley tree 
topology with G = 2 layers and coordination number z = 3 (one mother and k = 2 daughter 
branches). Dendritic branchlets are driven by independent Poisson stimuli (small arrows), (c) 
Each dendritic branchlet can be in one of three states: quiescent (0), active (1) or refractory (2). 
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Each site represents a dendritic branchlet, which we model with a three-state excitable 



element 



36| : Xi{t) G {0,1,2} denotes the state of site i at time t. If the branchlet is 



active (xj = 1), in the next time step it becomes refractory (xj = 2) with probability ps- 
Refractoriness is governed by p^, which is the probability with which sites become quiescent 
(xj = 0) again [see Fig. [U^c)]. Here we have used ps = I and p^ = 0.5. The propagation 
of dendritic spikes along the tree is assumed to be stochastic as well: each active daughter 
branch can independently excite its mother branch with probability px, contributing to 
what is referred to as forward propagation [i.e., from distal dendrites to the soma, see large 
descending arrow in Fig. [IJb)]. Backpropagating activity is also allowed in the model, 
with a mother branch independently exciting each of its quiescent daughter branches with 
probability (3px [large ascending arrow in Fig. [T](b)], where < /3 < 1. 

Dendrites are usually regarded as the "entry door" of information for the neuron, i.e., the 
dominant location where (incoming) synaptic contacts occur. Our aim then is to understand 
the response properties of this tree-like excitable medium. Incoming stimulus is modeled 
as a Poisson process: besides transmission from active neighbors (governed by px and 
each quiescent site can independently become active with probability ph = 1 — exp{—hAt) 
per time step [see Fig. [I^c)], where At = 1 ms is an arbitrary time step and h is referred 
to as the stimulus intensity. It reflects the average rate at which branchlets get excited. 



20|. With 



after the integration of postsynaptic potentials, both excitatory and inhibitory 
synchronous update, the model is therefore a cyclic probabilistic cellular automaton. 

A variant of the model accounts for the heterogeneous distribution of synaptic buttons 



along the proximal-distal axis in the dendritic tree. It consists of a 
h{g) = /ioc"^, with a controlling the nonlinearity of the dependence 
restrict ourselves to the simpler cases /3 = 1 and a = 0. 



ayer-dependent rate 



20|. We will mostly 



A. Simulations 

In the simulations, the activity F of the apical {g = 0) dendritic branchlet is determined 
by the average of its active state over a large time window (T = 10^ time steps and 5 
realizations). The response function F{h) is the fundamental input-output neuronal trans- 
formation in a rate-code scenario (i.e., assuming that the mean incoming stimulus rate and 
mean output rate carry most of the information the neuron has to transmit). 
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In a never ending matter of investigation, rate code has historically competed with tem- 



poral code, which is also supported by plenty of evidence 
tection 
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37|. Auditory coincidence de- 



as well as spatial localization properties of place and grid cells fundamentally 



depend on the precise spike time 



391]. Spike-timing-dependent plasticity, responsible for 



memory formation and learning, critically relies on small time differences (of order of tens 
of milliseconds) between presynaptic and postsynaptic neuronal spikes 40|. Moreover, zero- 
lag or near zero-lag synchronization, which are thought to play an active role in cognitive 
tasks 4l|, has been recently shown to be supported and controlled by neuronal circuits de- 



spite long connection delays 



42 



-|45|. Nevertheless, because of its robustness to the high level 



of stochasticity and trial-to-trial variability present in the brain 46|, rate code is probably 



more globally found 
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In this paper we implicitly assume that rate code holds. 



A typical response curve obtained from simulations with px = 0.7 and G = 10 is shown 
in Fig. [2] (symbols). It is a highly nonlinear saturating curve, with the remarkable property 
of compressing decades of stimulus intensity h into a single decade of apical response F. A 
simple measure of this signal compression property is the dynamic range A, defined as 



A = 10 log 



'•90 
^-10, 



(1) 



where = F ^{Fr^) is the stimulus value for which the response reaches x% of its maximum 
range: F^ = Fmin + ^{F^ax - i^mm), where F^^n = lim,,^o F{h) and F^a 



= lim/^^oo^(/i)■ 
As exemplified in Fig. Ej A amounts to the range of stimulus intensities (measured in dB) 
which can be appropriately coded by F, discarding stimuli which are either so weak as to be 
hidden by the self-sustained activity of the system [F < Fiq) or so strong that the response 
is, in practice, non-invertible owing to saturation {F > Fqq). 



Several features of this model have been explored previously 20j, like the dependence of 
A on model parameters, the double-sigmoid character of the response function, as well as the 
robustness of the results with respect to variants which increase biological plausibility. All 
these were based on simulations only. We now attempt to reproduce the results analytically. 




FIG. 2. (Color online) Response curve F{h) for simulations (symbols) and mean-field approxima- 
tions (lines; IS, 2S and black for EW) for px = 0.7 and G = 10. Horizontal and vertical arrows 
show the relevant parameters for calculating the dynamic range A (see Eq. [T|) . 

III. MEAN-FIELD APPROXIMATIONS 



A. Master equation 

The system can be formally described by a set of master equations. For the general 
case of arbitrary k, let Pf {x]y]u^^\v^^\ . . }j be the joint probability that at time t a site 
at generation g is in state y, its mother site at generation (7 — 1 is in state x, i (j) of its 
daughter branches at generation g + 1 are in state u (v) etc. 

Although the results in this paper are restricted to trees with k = 2, for completeness 
we write down the master equation for general k. The explicit derivation of the master 
equations for any layer is shown in Appendix |Al The equations for < g < G can be 
written as follows: 
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(l-p,)Pf(;l; 



P4i(; 2; ) = psPii; 1; ) + (1 - p,)P/(; 2; 
P/+i(;0;) = l-P4i(;l;)-P/+i(;2;), 



(2) 
(3) 
(4) 



where P{'{x;y;w^'^^) = P/(x;?/;) is a two-site joint probability and Pt{;y;) [also written 
Pt{y) for simplicity] is the probability of finding at time t a site at generation g in state y 
(regardless of its neighbors). 

Equations for the central {g = 0) and border [g = G) sites can be obtained from straight- 
forward modifications of Eq. |2l rendering 



A°+i(;i;) = A°(;0;)-(i-p.)5: 

4 = 

+ (1-P,)P0(;1;), 



P,^i(;l;) = P,^(;0;) + (l-p,) /3paP^'(1; 0; ) + (1 - P5)P,^(; 1; 



(5) 
(6) 



whereas Eqs. |3] and H] remain unchanged. Naturally, the full description of the dynamics 
would require higher-order terms (infinitely many in the limit G — t- oo), but Eqs |2] to H] 
suffice to yield the mean-field equations we address below. 

B. Single-site mean-field approximation 

The simplest method for truncating the master equations is the standard single-site (IS) 



49| . which results from discarding the infiuence of any neighbors in 

Piiy)- If this procedure is 



9+1/ 



(15) 



mean-field approximation 

the conditional probabilities: Pi{y\x) = Pi{;y;x)/P^ 

applied separately for each generation g, one obtains the factorization P/ (^x; y; u 
Pi~^{x)P^ {y)[Pi^^{u)Y[Pi^^{y)]^ , which reduces the original problem to a set of coupled 



(IS) 



equations for single-site probabilities: 



(IS) 



P/+i(l) VP/(0)A^(t) + (l-p5)P/(l) 



where 



A^'(t) 



Ph) 



i-(3p,pr\i)] \i-pxPr\i) 



(7) 



(8) 



is the probability of a quiescent site becoming excited due to either an external stimulus or 
propagation from at least one of its z = k + 1 neighbors (i.e., for < g < G). For g = and 
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g = G one has 



r 1 k-\-l 

A'{t) = l-{l-p,) l-p,Pl{l)] , (9) 



A^(t) = 1 - (1 -p,) [1 - (3p,P,^-\l)\ . (10) 

Note that this approximation retains some spatial information through its index g, whereby 
the generation-averaged activation -Pf (1) is coupled to Pf"^^(l) and Pf~^(l), rendering a 
2{G + l)-dimensional map as the reduced dynamics [note that the dimensionality of the 
probability vector is 3(G + 1), but normalization as in Eq. |4] reduces it to 2(G + 1)]. Although 
this facilitates the incorporation of finite-size effects (which are necessary for comparison 
with finite-G system simulations), we will see below that the results are not satisfactory. In 
fact, the results are essentially unchanged if we further collapse the different generations: 
Pt{x) = Pt{x), Vf? (which is the usual mean- field approximation, implying surface terms are 
to be neglected in the limit G — t- oo). The reasons for keeping a generation dependence will 
become clear when we propose a different approximation (see section UlI Pp . 

To compare our results with the case of interest for real dendrites, in the following we 
restrict ourselves to the binary tree, namely k = 2. Figure shows the results for the 
stationary value F = limj_^.oo -P°(l) in the absence of stimulus, i.e., the fixed point of the 
IS mean-field equations for /i = 0, as a function of the branchlet coupling px. The parameter 
values are G = 10, /3 = 1 and ps = 1 (deterministic spike duration). In the absence of 

(IS) 

stimulus, we see that the IS approximation predicts a phase transition at px = Pxc =1/3. 
As a consequence, the response curve F{h) for p > p\^ displays a plateau in the limit /;,—)■ 0, 
as shown in Fig. [3](b). The IS approximation yields results comparable to simulations only 
below Pxc 1 but performs rather poorly above the phase transition it predicts. However, 
given the deterministic spike duration (the only state in which a given site can excite its 
neighbors) and the absence of loops in the topology, a stable phase with stable self-sustained 



activity cannot exist 20L l24 1 . 



Figure [3](b) also shows the response curves as predicted by the simplified equations ob- 
tained from the G — )■ oo limit [i.e., by collapsing all layers, P/(x) = Pt{x), "ig]. Since they 
nearly coincide with the equations for G = 10 (which have a much higher dimensionality), 
it suffices to work with G — t- oo, which lends itself to analytical calculations. By expanding 
(around F ~ 0) the single equation resulting from Eqs. [3l HI [7] and [8] in their stationary 

10 




FIG. 3. (Color online) Mean-field approximations, (a) Firing rate of the IS, 2S and EW (black) 
approximations in the absence of stimulus {h = 0). (b) Family of response functions for (3 = 1 



201 ]) and curves are the 



and px = 0,0.2,0.4,...,!. Symbols represent simulations (as in Ref. 
IS mean-field approximation. Open (closed) symbols correspond to probabilistic (deterministic, 
Px = 1) neighbor coupling, and dotted (continuous) curves correspond to G = 10 (G oo). (c) 
Phase diagram under the IS approximation for different system sizes, (d) Same as (b) but for 2S 
approximation with G ^ oo. 

states, one obtains the value of critical value of px as predicted by the IS approximation for 
general k, ps and (3: 

(15) _ PS 

Still with Ph = [i.e., in the absence of stimulus {h = 0)], and recalling that At = 1 ms (i.e., 
rates F and h are expressed in kHz), the IS approximation yields the following behavior 
near criticality (i.e., for px ^ Pxf^)'- 

F(/^ = 0,e)^|6^ (12) 
11 



where /3 = 1 is a critical exponent, e = ^^_^ijs) i and C = ^^^^^2 ■^^^ + /? + ^^diE^.^ Since 



p 



Ac 



2 '^J 



in this case the order parameter corresponds to a density of activations and the system has 
no symmetry or conserved quantities, /3 corresponds to the mean-field exponent of systems 
belonging to the directed percolation (DP) universality class 49 1. 

The response function can also be obtained analytically for weak stimuli (for h <^ e, 
P/j <^ 1, thus ph — hAt = h). Below criticality {px < P^xf^), the response is linear: 

F{h,e)c^^. (13) 



As is usual in these cases, the linear response approximation breaks down at px = p^xc^ 24 1. 
For Px = Pxf^ , one obtains instead 

F(/.,e = 0)^ (^-j , (14) 

where 5/i = 2 is again a mean-field exponent corresponding to the response at criticality 49|. 

In Fig. ^c) we show in the plane {px, /?) the critical line given by Eq. [TTl as well as 
the line obtained by numerically iterating Eqs. [TlfTUl for finite G. It is interesting to note 
that the curves for G — t- 00 and G = 5, 10 split when /3 decreases. If one remembers that 
the simulated model has no active phase, the resulting phase diagram suggests that the IS 
solution can perform well for /3 ^ 0. Unfortunately, however, the limit /3 — )■ corresponds 
to the absence of backpropagating spikes, which in several cases of interest is far from a 



realistic assumption (backpropaga 



been observed experimentally 
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;ion of action potentials well into the dendritic tree has 



5l|). 



C. Two-site mean-field approximation 

The next natural step would be to consider the so-called pair or two-site (2S) mean-field 

(25) 



approximation 



49| . in which only nearest-neighbor correlations are kept: Pt{x\y,u,v) 



Pt{x\y). In that case, the dynamics of one-site probabilities end up depending also on two- 



site probabilities (2J]. Those, on their turn, depend on higher-order terms, but under the 2S 
truncation these can be approximately written in terms of one-site and two-site probabilities. 
The schematic representation of a general pair of neighbor sites {x and y), along with their 
corresponding neighbors (a, b and u, v), is depicted in Fig HI In the case of an infinite tree, 
and restraining oneself to the isotropic case /? = 1 , one can drop the generation index g and 

12 



employ the isotropy assumption Pt{x,y) = Pt{y,x) to write the general joint probability in 
the two-site approximation as 

P,{a; X- y, b; u, v) ^ [P,{x)P,{y)Y " ^^^^ 




FIG. 4. Two-site mean-field approximation schematic representation of a general pair [x and y) in 
a binary tree {k = 2). In order to describe the dynamics of x and y, each of their neighbors must 
be taken into account. According to Eq. [T5l the joint probability of the labeled sites is rewritten 
in terms of two-site probabilities. 

In this simplified scenario, the collective dynamics is reduced to that of a probability 
vector containing two-site probabilities (from which single-site probabilities can be obtained, 
please refer to Appendix |B]). Taking all normalizations into account, the dimensionality of 
this vector can be reduced to 5. As can be seen in Appendix [Bl however, this simple 
refinement in the mean-field approximation already leads to very cumbersome equations. 

As shown in Figs. El El^a), and[3t^d), the gain in the quality of the approximation falls far 
short of the increase in the complexity of the calculations. In fact, the IS and 2S approx- 
imations yield qualitatively similar results, capturing the essential features of the system 
behavior only for p\ smaller than some critical value p\c. For px > pxc, both approximations 
predict a phase transition to self-sustained activity, with pxc = 1/3 for IS and pxc = 1/2 for 
2S (in the case /3 = 1 = ps). These predictions are incorrect: when simulating the model 
without external driving (/i = 0), in a few time steps [(9(6*)] the system goes to an absorbing 



state 491], from which it cannot scape in the absence of further stimulation. 

One can interpret the results of the approximations as follows. At the IS approximation 
level, a quiescent site will typically be activated by any of its three spiking neighbors at 
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the phase transition, hence pxc = 1/3. The refinement of the 2S approximation consists 
in keeping track of the excitable wave propagation from one neighbor, leaving two other 
neighbors (wrongly assumed to be uncorrelated) available for activity propagation, hence 
Pxc = 1/2. 

One could, in principle, attempt to solve this problem by increasing the order of the 
cluster approximation (keeping, e.g., 3- and 4-site terms). However, the resulting equations 
are so complicated that their usefulness would be disputable, especially for applications in 
Neuroscience. It is unclear how more sophisticated mean- field approaches (such as, e.g.. 



non-equilibrium cavity methods 



52|-|54|) would perform in this system. In principle, they 



seem particularly appealing to deal with the case ps < 1, when a phase transition to an 
active state is allowed to occur (and whose universality class is expected to coincide with 



55 



56|). Attempts in this direction are promising and 



that of the contact process on trees 
would be welcome. 

In the following section, we propose an alternative approximation scheme which circum- 
vents the difficulties of the regime px I and at the same time takes into account finite-size 
effects. 



D. Excitable- wave mean-field approximation 

The difficulties of the IS and 2S approximations with the strong-coupling regime are not 
surprising. Note that the limit of deterministic propagation (approached in our model as 
Px ^ 1) of deterministic excitations {ps = 1) is hardly handled by continuous-time Markov 
processes on the lattice. To the best of our knowledge, a successful attempt to analyti- 
cally determine the scaling of the response function to a Poisson stimulus of a hypercubic 
deterministic excitable lattice was published only recently js^] (and later confirmed in bio- 



physically more detailed models [29|). While these scaling arguments have not yet been 
adapted to the Cayley tree, the collective response resulting from the interplay between 
the propagation and annihilation of quasi-deterministic excitable waves remains an open 
and important problem. In the following, we restricy ourselves to the case ps = 1, i.e., 
deterministic spike duration. 

As discussed above, the IS and 2S approximations give poor results essentially because 
they fail to keep track of where the activity reaching a given site comes from. We therefore 
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propose here an excitable-wave (EW) mean-field approximation which attempts to address 
precisely this point. 

The rationale is simple: in an excitable tree, activity can always be decomposed in 
forward- and backward-propagating excitable waves. Formally, this is implemented as fol- 
lows. We separate (for g > 0) the active state (1) into three different active states: lA, IB, 
and IC, as represented in Fig. [5t^a). (lA) stands for the probability that activation (at 
layer g and time t) was due to the input received from an external source (controlled by ph). 
The density of elements in lA can excite quiescent neighbors at both the previous and the 
next layers. (15) corresponds to the density of elements in layer g which were quiescent 
at time t — 1 and received input from the next layer {g + 1) (i.e., a forward propagation). 
The density of elements in IB can excite solely quiescent neighbors at the previous layer. 
Finally, P/(1C) accounts for the activity coming from the previous layer (i.e., backpropa- 
gation). The density of elements in IC can excite solely quiescent neighbors at the next 
layer. For lack of a better name, we refer to these different virtual states as excitation 
components. Figure O^b) represents the activity flux in the dendritic tree as projected by 
the EW mean-field approximation. The absence of loops guarantees the suppression of the 
spurious non-equilibrium phase transition predicted by the traditional cluster expansions. 

Following these ideas, one can write the equations for the g > layers as 

P^IA) = Pf (O)A^ , (16) 
P/+i(lP) = Pf(0)(l-A^)Af,(t), (17) 
PU^C) = Pf (0)(1 - A^) [1 - AUt)] m) , (18) 

where, in analogy with Eq. [HI the excitation probabilities are now given by 

= Ph, (19) 

A%it) = i-{i-p, [pr\iA) + pr\iB)]Y , (20) 

AUt) = (3px \Pr\lA) + Pr\lC)] . (21) 



Equations E] and II remain unchanged, with P/'(l) = Pf (lA) + P/(1P) + Pf (IC). The dy- 
namics of the most distal layer g = G is obtained by fixing A^(t) = 0. The apical {g = 0) 
element has a simpler dynamics since it does not receive backpropagating waves, so its activ- 
ity is governed by Eq. U\ with ^ = and A%t) = 1 - (1 -p^) {1 - px [P/(1A) + P/(1P)]}^^^ 

15 




FIG. 5. (Color online) Schematic representation of the excitable-wave mean-field approximation, 
(a) Dynamics of each layer {g > 0) in the excitable-wave (EW) mean-field approximation (see 
text for details). There are 3 different active states: lA represents activity coming from an ex- 
ternal input; IB is reached due to forward activity from the next layer, whereas IC is excited by 
backpropagating activity from the previous layer, (b) Schematic EW mean-field approximation 
dynamics in a tree with G layers. Note that there are no loops in the activity flux. 



instead of Eq. |8l Taking into account the normalization conditions, the dimensionality of 
the map resulting from the EW approximation is 4{G — 1) + 5. 

It is important to notice that, while Eqs. [T91I2T] are relatively straightforward, there is 
a degree of arbitrariness in the choice of Eqs. [T6lfT8l As written, they prescribe an ad hoc 
priority order for the recruitment of the excitation components of the EW equations: first 
by synaptic stimuli (Eq. [T6|) . then by forward propagating waves (Eq. [T7|) . and finally by 
backpropagating waves (Eq. [T5]) . This choice seems to be appropriate in the regime of 
weak external driving, insofar as the order coincides with that of the events observed in 
the experiments: forward dendritic spikes, a somatic spike, then backpropagating dendritic 



spikes 5l| • Appendix [C] compares the response functions for different priority orders to 



emphasize the robustness of the approximation with respect to that. 
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FIG. 6. (Color online) Response functions: simulations compared to the EW approximation (both 
with 13 = 1) and experimental data, (a) Family of response functions for G = 10 and px = 



201 1) and solid curves are the EW mean- 



0, 0.2, 0.4, . . . , 1. Symbols are the simulations (as in Ref. 
field approximation, (b) Family of response functions for px = 0.7 and different tree sizes: G = 
5, 10, 15, 20. (c) Simulations and EW response function for external stimuli spatially distributed 
as h{g) = /loe''^: from right to left, a = 0, 0.1, 0.3, . . . , 0.9 (px = 0.8 and G = 10). Horizontal lines 
are plotted for the estimation of the dynamic range, (d) Experimental result from mouse retinal 



ganglion cells shows double sigmoid response curves Jc 
(measured in rhodopsins/isomerizations/rod/second 



osed symbols) function of stimulus / 



58(1). They can be reasonably well fit by both 



simulations (open symbols, px = 0.58 and h = 0.37 I) and the EW response function {px = 0.59 
and h = 0.4 I) with tree size G = 15. 



Though noncontrolled, the EW mean-field approximation does provide excellent agree- 
ment with simulations. The results for G = 10 can be seen in Fig.[6]^a), which shows a family 
of response curves F{h) for varying coupling px- One observes that the EW mean- field re- 
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suits (lines) follow the simulation results (symbols) very closely up to px — 0.8, reproducing 



even the double-sigmoidal behavior of the curves |20l. |58N60| . For larger values oi px, agree- 
ment is restricted to very small or very large values of h (for intermediate values of h, note 
that F( p\) is nonmonotonous, a rather counterintuitive phenomenon called "screening reso- 



nance" 



20|). Most importantly, however, the EW equations eliminate the phase transition 
wrongly predicted by the traditional mean-field approximations. 

In a real neuron, the number of layers is finite [0(10) or less] and it would be extremely 
interesting to have an analytical approximation which managed to take finite-size effects 
into account. As it turns out, the mean- field approximation we propose can do precisely 
that, since it couples densities at different layers (so G again controls the dimensionality of 
the mean- field map). Figure [^b) compares simulations (symbols) with the stationary state 
of the EW mean-field equations (lines) for different system sizes. Note that the agreement 
is excellent from G = 5 up to G = 20, for the whole range of h values. 

The EW mean-field approximation is also very robust against previously proposed vari- 
ants of the model. For instance, in several neurons the distribution of synaptic inputs 
along the dendritic tree is nonuniform, increasing with the distance from the soma. A one- 
parameter variant which incorporates this nonuniformity consists in a layer-dependent rate 



h{g) = hoe""^ [20|, as described in section [TTl Figure EJ^c) depicts a good agreement between 
simulations and the EW approximation for a range of a values. 

Also shown in Fig. M^d) is a comparison among experimental results from retinal ganglion 



cells to varying light intensity [58| (closed symbols), simulations (open symbols), and EW 
mean- field approximation (lines), which agree reasonably well. Therefore the approximation 
we propose can, in principle, be useful for fitting experimental data and reverse-engineering 
parameter values from data-based response functions at a relatively small computational 
cost (say, compared to simulations). In this particular example, it is important to emphasize 
that the experimental response curves are, in principle, influenced by other retinal elements. 
Given the very simple nature of our model, it is hard to pinpoint which part of the retinal 



circuit our Cayley tree would represent. Following Shepherd |20|, |61[, however, we suggest 
that the ganglionar dendritic arbor plus the retinal cells connected to it by gap junctions 
(electrical synapses) can be viewed as an extended active tree similar to the one studied 
here, with a large effective G. 

The dynamic range A is one of the features of the response function which has received 
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attention in the literature in recent years 
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33 



35| . Here it serves the purpose of summariz- 



ing the quahty of the EW mean-field approximation in comparison with model simulations. 
In Fig. [7] we plot A as a function of px for several system sizes G. Both the IS and 2S 
mean-field approximations predict a non-equilibrium phase transition in the model, where 
a peak of the dynamic range therefore occurs 25|]. Both approximations perform badly es- 



pecially in the high-coupling regime. The EW approximation correctly predicts the overall 
behavior of the A{px) curves, for all system sizes we have been able to simulate. Finally, 
the inset of Fig. [7] shows a second variant of the model in which the parameter /3, which 
controls the probability of a spike backpropagating, is free to change. Once more, the EW 
approximation manages to reproduce the A{px) curves obtained from simulations for the 
full range of (3 values. 








0.2 0.4 „ 0.6 0.8 



FIG. 7. (Color online) Dynamic range as a function of the coupling parameter for IS (G = 10), 
2S (infinite tree), EW (black lines) mean-field approximations compared to simulations (symbols, 



as in Ref. [20|) for different trees sizes G. The inset shows the dynamic range of dendritic trees 
subjected to an asymmetrical activity propagation probability controlled by /3 (symbols stand for 
the simulations and the curves for the EW approximation for G = 10). 
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IV. CONCLUSIONS 



The need for a theoretical framework to deal with active dendrites has been largely recog- 
nized. However, the plethora of physiological details which are usually taken into account to 
explain local phenomena renders the problem of understanding the dynamics of the tree as a 
whole analytically untreatable. What we have proposed is the use of a minimalist statistical 
physics model in which our ignorance about several physiological parameters is thrown into 
a single parameter px- The model has provided several insights and predictions, most of 
them yet to be tested experimentally 20(]. Here we have shown that the model is amenable 
to analytical treatment as well. 

We have compared different mean-field solutions to the model, and shown that standard 
cluster approximations (IS and 2S) yield poor results. They incorrectly predict phase transi- 
tions which are not allowed in the model, thereby failing to reproduce the response functions 
precisely in the low-stimulus and highly nonlinear regime (where most of the controversies 



are bound to arise |62| ) . 

To overcome this scenario we developed an excitable-wave mean-field approximation 
which takes into account the direction in which the activity is propagating through the 
different layers of the tree. Though ad hoc, the approximation reproduces simulation results 
with very reasonable accuracy, for a wide range of parameters and two biologically rele- 
vant variants of the model. We hope that our EW mean-field approximation may therefore 



contribute to the theoretical foundation of dendritic computation (63 1. 

It is important to recall that the theory attempts to address a model in a regime which is 
expected to be close to that of a neuron in vivo: dendritic spikes are generated at random, 
may or may not propagate along the dendrites, and annihilate each other upon collision. 
Our model allows one to formulate theoretical predictions for in vivo experiments in sensory 



system neurons, e.g., ganglion cells (58 



^olfactory mitral cells or their insect counterparts. 



i.e., antennal lobe projection neurons (62|: a) if one removes part of the dendritic tree and/or 
b) blocks the ionic channels responsible for dendritic excitability, one should see a decrease 
in the neuronal dynamic range. To the best of our knowledge, these experiments have not 
been done yet. 

Another issue upon which our results could have a bearing is the so-called gain control 
modulation. In the neuroscience literature, the term refers to the neuronal capacity to change 
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the slope of the input-output response function IgsJ . This property has been reported in 



visual 



66 



-l69j. somatosensory 



ion J65 
.ry Inl 



70[ and auditory [71[ systems. Several possible mechanisms 



have already been proposed to explain gain control, based on synaptic depression 
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background synaptic input 



73 



74( 1 . noise 75 L shunting inhibition 76|, |77| and excitatory 



(NMDA 78|) as well as inhibitory (GABA^ 79|) ionotropic receptor dynamics. 

All these mechanisms are intercellular, in the sense that they rely on the influence of 
factors external to the neuron. Our model, on the other hand, shows gain control in its de- 
pendence on the coupling parameter px, which controls the propagation of dendritic spikes 
within the tree. It is therefore an intracellular mechanism which offers an additional expla- 
nation for this ubiquitous phenomenon. 

The physics of complex systems is becoming more and more embracing, shedding light 
in different areas, including neuroscience. Particularly at the cellular and subcellular levels, 
we foresee the merging of the two fields, dendritic computation and statistical physics, 
as a promising avenue. The maturity of the latter could illuminate several frontiers in 
neuroscience. 
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Appendix A: Derivation of the general master equation 

Let us illustrate how the master equation is obtained by starting with the simplest possible 
case, namely, the latest layer of the Gayley tree. Sites at the surface connect to a single site 
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(their mother branchlet), so the probabihty of their being excited at time t + 1 is 

P,%,i; 1; ) = [1 - (1 - - Pp.m'^il; 0;)+phY: (/; 0; ) + (1 - ps)P['i; 1; ) • (Al) 

1^1 

Each term has a straightforward interpretation: the first term corresponds to the probabihty 
that the surface site is quiescent (state 0), its neighbor is active (state 1), and excitation 
gets to the surface via an external stimulus (ph) and/or backpropagating transmission (Ppx)] 
the second term corresponds to the excitation of the surface site via an external stimulus, 
provided that it is quiescent (0) and its neighbor is in any state other than active (1); the 
third term corresponds to the probability that the surface site was in state 1 and did not 
move to state 2 (a transition controlled by ps, see Fig. [It). 

The next easiest case to consider is that of the root {g = 0) site. It connects to + 1 
daughter branchlets, and can be excited by any number of them. Contrary to the g = G 
surface sites, the root site only receives forward propagating activity (hence /3 plays no role). 
Analogously to Eq. lAlt the equation for Pt^i is given by 



pO /.I. 



'k + V 

,k + l. 



l-(l-p,)(l-p,)'=+i]P,°(;0; 1(^+1)) 



.k-1. 



+ 



'k + V 



+ [1 - (1 - p,)(l - PxD E P'i; 0; Jl, J2, . . .,Jk+l-m) I 

jl :i2vijfc+l-m^l 
+ ■■■ 

+ [1-{1-Ph){1-Px)] E P,°(;0;l,ji,...,j.)f^t^) 
+Ph E Pt°(;0;ji,...,jfc+i) 

+ (l-p,)P,°(;l;). (A2) 

The terms of the kind [1 — (1 — ph){l — P\)"^] account for the excitation of the root site 
via an external stimulus and/or via transmission from m of its active daughter branch- 
lets, regardless of the state of its k + 1 — m non-active neighbors (hence the sum over 
Ji) J2! • • • ijk+i-m 7^ !)• Each term is weighted by the number of combinations of m 

active sites out of A; + 1. The sum of these terms (from m = 1 to m = k + 1) therefore plays 
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a role equivalent to that of the first term in eq. lAll The two last terms are analogous to 
those of eq. lAll 

Finally, we come to the equation for a general site with 1 < g < G — 1, which can be 
excited by both its mother as well as its daughter branchlets. The equation for Pt+i thus 
generalizes the terms of the preceding equations: 



+ 
+ 



' k ' 
.k~l. 



+ [!-(!- p,){i - - p,)] Pt{i; 0; 1, ji, . . . , j.-i) C^] 

ji,...,jk-i^i VJ 

7l,...,.7fc^l V^/ 



+ 

+ 
+ 



1-(1-P,)(1-Pa)'=1EA'(^;0;1('=)) 



.k, 



' k ^ 
.k-l. 



[1 - (1 - pn){i - px)] Y Pti^; 0; 1, ji, • • • , j.-i) i\ 



+Ph Y Pf(^;0;ji,...,jfc) 



(A3) 



Equations I AltlASl can be drastically simplified 2^- Taking into account the normalization 
condition 

P/(a; b;ji, . . . , = ^Pf (a; b;ji, . . , (A4) 

jc 

the sums in eqs. IAltlA3l can be reduced. For instance, 
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Pf(£;0;l('=-^)) -P/(£;0;l('=)) . 



(A5) 



Iterating this procedure and rearranging terms, one finally arrives at 



P,%(;l;)=P°(;0;lW)-(l-j,.) 



-1W(;0;1«) 



(A6) 
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P,%(;l;) = P/(;0;lW)-(l-p.) 



E(k(^)(-iW(;0;1«) 



+ (l-p,)P/(;l;), (AT) 



Pt%ii: 1; ) = 0; 1^°)) - (i - ph) [p?{: O; 1^°)) - PpxP?{i-, O; 

+ (l-p5)Pf(;l;). (A8) 
Recalling that P/(a;; = P/(x; ), we recover eqs. [2], |5]and[6l 

Appendix B: Two-site mean-field equations 

For an infinite dendritic tree, the complete set of equations under the 2S approximation 
is given by: 



(25) 

Pi+i(0; 0) {(1 - p,)'Pi(0; Q)[pxPt{^; 1) - Pi(0)]^ 

+2p^{l - pH)PtiO; 2)PtiQy[pxPtiO; 1) - PtiO)]^ 
+p5P,(2;2)P,(0)n^; (Bl) 



Pi+i(0; 1) -(Pi(0; l)Pi(0)2{-(l - - ph)il - Px)[pxPtiO; 1) - PtiO)]'} 
+ (1 - p,)P<(0; 0){(1 - Ph)[pxPt{0; 1) - Pt(0)]^ 

-Ptmpxm-A) - pmn 

+p,PtiO; 2){(1 - p,)Pi(0)2bAP<(0; 1) - PtiO)]' - ^^(O)'} 

-(1 -p,)p^P,(l;2)P,(0r)-^ ; (B2) 

Pt+i(0; 2) {(1 - - p,)Pi(0; 2)[pAPt(0; 1) - Pt{0)f 

+ps{l-ph){l-px)Pt{0; l)[pxPt{0; 1) - Pt(0)]' 
+p^(l-p^)Pi(2;2)Pi(0)2 

+p,p^P,(l;2)Pi(0)n-^; (B3) 
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-2P,{of\pxm, 1) - pmr + pm'} 

-2Pi(0; l)Pt{Of{il-ps){l-Ph){l-px)[pxPtiO; 1) - PA^''' 



-(1 - ps)Pt{or} + (1 - psfPtii; mior) ; (b4) 



Pt+i{l; 2) -(P,(0; 2){(1 - p,){l - PH)[pxPt{0; 1) - Pt{0)? - (1 - P7)^t(0)'} 
+PsPt{0; 1){(1 - - Px)[pxPtiO; 1) - Pi(0)]2 - P,(0)2} 
-p5(l-P5)Pi(l;l)Pt(0)' 



:i -p^)P,(l;2)P,(0)2)-^ ; (B5) 



(2S) 

Pi+i(2;2) ^ p2p^(l;l) + 2p5(l-|)^)Pi(l;2) + (l-p^)2p,(2;2) . (B6) 

In order to solve it numerically, by normalization we make use of: 

Pi(0) = Pi(0; 0) + Pi(0; 1) + Pt(0; 2) . (B7) 

To compare the results with the simulations and experiments we use F = limf__j.oo Pt(l), 
where: 

Pi(l) = Pt{0; 1) + Pt{l; 1) + Pt{l; 2) . (B8) 
Finally, by completeness, the last single-site probability can be obtained by: 

Pt{2) = Pi(0; 2) + Pt{l; 2) + Pi(2; 2) . (B9) 

Appendix C: Robustness with respect to the priority order of the excitable com- 
ponents 



The priority order used in Eqs. [T6lfT8] was "ABC", i.e., first lA, followed by IB and then 
IC. The virtual state IB accounts for the forward-propagating excitable- wave flux whereas 
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state IC accounts for the backward- propagating excitable- wave flux. In order to compare 
all the different combinations, Fig. |8] displays families of response functions. 

Intuitively, the neuronal flring rate F increases as the forward-propagating excitable- 
wave flux grows. Summarizing the results, switching the order of IB and IC (lower panels 
in Fig. [H]), we reduce the forward-propagating excitable- wave flux, and consequently, the 
response functions corresponding to large p\ values present a lower flring rate. Moreover, 
the result is virtually the same irrespective of the order in which lA appears (compare panels 
horizontally in Fig. [8]). 




10'^ 1 10^ lO"* 10'^ 1 10- 10"* 1 10- 10"* 

h(s-') h(s-') li(s-') 



FIG. 8. Effects of the priority order of the excitable components on the response functions: simu- 
lations compared to the EW approximation and experimental data. Family of response functions 
for G = 10 and px = 0,0.2,0.4,... ,1. Top panels display combinations of IB prior to IC, and 
bottom panels the converse. Symbols are the simulations and solid curves are the EW mean-field 
approximation . 

The approximation is robust with respect to the order chosen, and only minor differences 
can be found for strong coupling ~ 1) and the intermediate amount of external driving 
input. Changes in the order of the components modify the response functions quantitatively. 
However, qualitatively, the response functions are very alike. Most of the differences occur 
at the high-coupling regime {p\ ~ 1), but do not affect the localization of hio and /igo (which 
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implies that the dynamic range remains unchanged). 
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